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Molecular dynamics (MD) methods are opening new opportunities for simulating the 
fundamental processes of material behavior at the atomistic level. However, increasing the 
size of the MD domain quickly presents intractable computational demands. A robust 
approach to surmount this computational limitation has been to unite continuum modeling 
procedures such as the finite element method (FEM) with MD analyses thereby reducing the 
region of atomic scale refinement. The challenging problem is to seamlessly connect the two 
inherently different simulation techniques at their interface. In the present work, a new 
approach to MD-FEM coupling is developed based on a restatement of the typical boundary 
value problem used to define a coupled domain. The method uses statistical averaging of the 
atomistic MD domain to provide displacement interface boundary conditions to the 
surrounding continuum FEM region, which, in return, generates interface reaction forces 
applied as piecewise constant traction boundary conditions to the MD domain. The two 
systems are computationally disconnected and communicate only through a continuous 
update of their boundary conditions. With the use of statistical averages of the atomistic 
quantities to couple the two computational schemes, the developed approach is referred to as 
an embedded statistical coupling method (ESCM) as opposed to a direct coupling method 
where interface atoms and FEM nodes are individually related. The methodology is 
inherently applicable to three-dimensional domains, avoids discretization of the continuum 
model down to atomic scales, and permits arbitrary temperatures to be applied. 


I. Introduction 

The emerging field of nanomechanics is providing a new focus in the study of the mechanics of materials: That of 
simulating fundamental atomic mechanisms involved in the initiation and evolution of damage. These simulations 
are commonly based on either quantum mechanics (ab-initio, tight-binding (TB), density-functional theory (DFT)) 
methods or on classical molecular dynamics (MD) and molecular statics (MS) methods. Predictions of material 
behavior at nanometer length scales promise the development of physics-based 'bottom-up' multiscale analyses that 
can aid in understanding the evolution of failure mechanisms across length scales. However, modeling atomistic 
processes quickly becomes computationally intractable as the system size increases. With current computer 
technology, the computational demands of modeling suitable domain sizes (on the order of hundreds of atoms for 
quantum mechanics based methods, and potentially billions of atoms for classical mechanics based methods) and 
integrating the governing equations of state over sufficiently long time intervals quickly reaches an upper bound for 
practical analyses. In contrast, continuum mechanics methods such as the finite element method (FEM) provide an 
economical numerical representation of material behavior at length scales in which continuum assumptions apply. 
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The concept of bridging length scales by coupling different computational paradigms is particularly attractive as 
a highly efficient means of reducing the computational cost of the simulations in cases that require modeling of 
relatively large material domains to capture the complete deformation field, but where atomic and subatomic 
refinement is needed only in very localized regions (e.g. near a crack tip or dislocation core). Such computational 
issues arise in modeling crack nucleation and propagation, and in modeling dislocation formation and interaction. 
By using coupled models, the size limitations of the atomistic simulation can be avoided by embedding an inner 
atomistic domain where complex dynamic processes and large deformation gradients exist within an outer domain 
where the deformation gradients are small so that a continuum representation of the material becomes appropriate. 

Over the past decade, various methods that couple material representations at different levels of refinement have 
been developed and offer significant computational advantages compared with full MD simulations for predicting 
deformation and fracture processes 1 " 11 . In all of these methods, the most challenging problem is the computational 
connection of the two different material representations at an imaginary interface where the continuity of the 
material properties and kinematics must be preserved during simulation. For example, the coupling between the 
subatomic quantum mechanics description and the classical atomistic mechanics representation requires a transition 
from a highly interdependent representation of the state of the whole system, which cannot be partitioned into 
individual atomic states, to the classical representation of individual atoms interacting through Newtonian forces. A 
way to resolve this challenge was provided by the work by Abraham and co-workers 1, 2 who developed a procedure 
that combines semi-empirical TB calculations, that account for the electronic structure as parameterized from first- 
principles of quantum mechanics, with a MD simulation based on Newtonian mechanics. The procedure was used 
to simulate the intricate mechanisms of atomic debonding at the tip of a propagating crack in silicon. 

No less challenging is the coupling between atomistic and continuum material representations. In this case, 
continuity of material properties must be maintained while transitioning from individual atoms interacting through 
nonlocal forces to the local stress-strain field formalism of continuum mechanics. Because of the potentially 
enormous computational benefit that might be achieved by replacing portions of the discrete atomistic domain with 
a continuum representation, allowing for orders of magnitude increase in the system size without a substantial 
increase of the computational cost, the atomistic - continuum coupling is attracting increased attention. Various 
specialized coupling methods 3 ' 5 have been developed to address different problems (e.g. crack propagation, 
dislocation evolution) in a computationally efficient manner. 

In one formulation, named the coupling of length scales method (CLS) 2 , the nodes in a finite element model 
(FEM) representing the continuum region are directly connected to the atoms of an atomistic region forming an 
interface of “pad” atoms. In another formulation, called the quasicontinuum (QC) method, reviewed by Miller and 
Tadmor 6 , the regions to be represented as either atomistic or continuum are determined by evaluating the magnitude 
of evolving local deformations. The representations of different regions are performed using representative atoms or 
“repatoms”. In the QC formalism, “nonlocal repatoms” are used to represent “real” atoms to form atomistic regions 
treated by MS/MD methods while “local repatoms” are used to define continuum domains using FEM procedures. 
Another representative coupling approach is the bridging method of Xiao and Belytschko 7 which is based on an 
overlay approach in which MD and FEM representations are superposed. The coupled atomistic/dislocation 
dynamics (CADD) method of Shilkrot et al . 8 is specifically designed to simulate, identify and pass dislocations 
between atomistic and continuum domains. For crack problems, the early efforts of Gumbsch and Beltz 0 led to the 
development of the finite element-atomistic (FEAt) coupling procedure that combined an embedded MD region with 
a finite element domain. A generalized formulation of the conventional FEM, which allows FEM nodes to be 
considered as coarse-grained MD “atoms”, led to another computational scheme for atomistic-continuum coupling 
called coarse grained molecular dynamics (CGMD). A detailed representation of CGMD is given by Rudd and 
Broughton 10 1 1 . Other variations and modifications of the above referenced approaches also exist. 

A common feature of all these approaches 1 " 11 for coupling atomistic with continuum representation is the 
refinement of the FEM mesh down to atomic scale and make a transition from FEM mesh nodes to discrete atoms at 
the interface. The approaches of relating the kinematics of atoms and FEM nodes in a one-to-one manner will be 
referred to in this paper as direct coupling (DC) approaches. 

While DC approaches are straightforward, the fundamental difficulty in their development lies in the inherent 
differences between the atomistic and continuum computational models. The physical state of the atomistic region 
is described through nonlocal interatomic forces between discrete atoms of given position and momentum, while the 
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physical state of the continuum region is described 
through continuous stress-strain fields that reflect local 
statistical averages of atomic interactions at larger 
length and time scales. Using a continuum 
representation at atomic length scales creates potential 
difficulties, such as the appearance of “ghost” forces 3 * " 5 
or the imposition of nonphysical constraints on the 
motion of atoms at the interface that are equated to the 
displacements of the continuum mesh rather than to the 
dynamical effects of interatomic forces. In general, the 
formal connection between continuum and discrete 
quantities can only be achieved through an adequate 
statistical averaging over scales where the discreteness 
of the atomic structure can be successfully 
approximated as a continuum. 


In this paper, an alternative approach to the 
DC approaches is proposed to construct a coupled MD- 
FEM system. The approach is based on solving the 
special boundary value problem (BVP) at the MD-FEM 
interface for a MD system embedded in a FEM domain 
as depicted in Figure 1. 

The method uses statistical averaging over both time and volume of atomistic subdomains at the MD-FEM 
interface to provide nodal displacement boundary conditions to the continuum FEM domain, which, in turn, 
generates interface reaction forces that are applied as constant traction boundary conditions 12 to the atoms within the 
localized MD subdomain. Thus, this approach may be described as a local-nonlocal BVP because it relates local 
continuum nodal quantities with nonlocal statistical averages of atomistic quantities over selected atomic 
subdomains. An iterative procedure between the MD statistical displacements and the FEM reaction forces ensures 
continuity at the interface. In this way, the problem of redefining continuum variables at the atomic scale is 
avoided, and the developed interface approach links different time and length scales between the MD and FEM 
domains. Typically, one finite element at the interface encompasses a region of several hundred to several thousand 
atoms. At this scale, the discreteness of the atomic structure is homogenized enough so that the FEM domain 
responds to the atomistic domain as an extension of the continuum. At the same time, the piecewise constant traction 
boundary conditions of the MD domain ensure that the elastic field from the FEM domain is correctly transferred to 
the atomistic region. With the use of statistical averages to couple the two computational schemes, the developed 
approach is referred to as a statistical coupling (SC) approach as opposed to the DC approaches. Based on the SC 
approach, the developed MD-FEM coupling method is referred to as the embedded statistical coupling method 
(ESCM). 

This paper will detail the ESCM approach to coupling MD and FEM computational domains for the case of 
systems that reach thermodynamic equilibrium or evolve quasistatically. While there is no principle difficulty in 
implementing this approach for non-equilibrium systems, it is beneficial to consider the case of equilibrium 
simulations first to illustrate the methodology. 

Section II. A describes the structure of the coupled MD-FEM model. Sections II. B and II. C discuss elements of 
the coupling interface involving the Interface MD Region and Surface MD Region, respectively. Section II. D 
presents an overview of the overall iterative MD-FEM coupling methodology. Section III presents several 
validation studies to substantiate the accuracy of the developed coupling methodology. Finally, Section IV will 
present concluding remarks on the overall effectiveness of the ESCM. 

II. The ESCM Model 

The ESCM model consists of four regions: the inner MD region, the interface MD region, the surface MD region 
and the FEM mesh as shown in Figure 2. The inner MD region consists entirely of an atomistic representation of 
the material. The interface MD/FE region consists of the interface MD region and the interface FEM nodes and 
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performs the primary coupling of MD and FEM 
domains via the ESCM statistical interfacing 
procedures. The surface MD region extends to the 
outer boundaries of the atomistic simulation and 
acts to approximate the internal forces of an 
infinite MD domain. Together, the inner MD 
region, the interface MD region and the surface 
MD region comprise the MD domain. Likewise, 
the FEM domain consists entirely of a continuum 
representation of the material and extends from the 
interface FEM nodes to the outer boundaries of the 
simulation. 


A. MD and FEM model components 

The inner MD region is used to model 
material phenomena at the atomistic level and is 
simulated under the full influence of the FEM 
reaction forces applied as traction boundary 
conditions to the interface region. This inner MD 
region should be large enough to sustain any of the 
types of processes (mostly plastic deformation 
events such as dislocation or void nucleation, 
crack propagation, etc.) that are required by the 
simulation. 

It is important to emphasize that the partitioning of the MD system into an inner region, an interface region and 
a surface region is not a physical separation of the system. An atom assigned to a particular location freely interacts 
with atoms in its interaction neighborhood that may reside in a different region. Thus, the overall simulation is 
performed using any conventional MD technique and the whole MD configuration is treated as an independent 
integrated MD system without any imposition of direct kinematic constraints. The only difference between the 
three MD regions is that, while the atoms in the inner MD region are subject only to their interatomic forces, the 
interface and surface MD regions are used to facilitate the application of external forces involved in ESCM coupling 
that act in addition to the interatomic forces. The modification of the MD computations relates only to the 
application of external forces and will be discussed in Sections II. B and II.C where details of the interface and 
surface MD regions are presented. 

The FEM domain is used to replace an atomistic representation with a continuum model in parts of the domain 
where the deformation gradients are small and atomic level resolution is not necessary. Using a FEM model 
permits a large reduction in the computational cost of performing MD simulations by limiting the atomic resolution 
to nonlinear processes in the inner MD region. The current application uses the FEM domain to simulate an 
extended material model such that the elastic deformation and load transfer due to applied far-field boundary 
conditions are accurately transferred to the inner MD region. The continuum field is currently assumed to be static 
with linear elastic material properties but other applications of ESCM might require the incorporation of nonlinear 
material behavior such as plasticity or general dynamic response such that nonlinear processes generated in the inner 
MD domain could be propagated into the continuum. 

B. Interface MD region 

The main role of the interface MD region (Figure 2) is to provide the connection between the MD and the FEM 
domains. The interface MD region is formed by atoms that surround those FEM nodes, called interface FEM nodes, 
that overlap the MD domain. The atoms that surround a given interface node form a cell in the interface MD region, 
called an interface volume cell (IVC). The IVCs are established to communicate averaged MD displacements 
computed at their mass center to the FEM nodes, and to provide a set of atoms over which FEM reaction forces are 
distributed to the MD system during simulation. The IVCs need not coincide in size and shape with the finite 
element to which the FEM node belongs. In the model described in this paper, the IVCs are formed through a 
Voronoi-type construction 13 by selecting those atoms with a common closest FEM node. Accordingly, the average 
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size of the IVCs becomes equal to the average distance between the interface FEM nodes, /, which also defines the 
in-plane thickness of the interface MD region (Figure 2). 

During subsequent coupled MD-FEM simulation, the volume of the IVC is the integration volume for the 
statistical displacement vector 5 j , that is prescribed to the interface FEM node. The displacement vector for each 

interface FEM node is calculated as a time-average of the center of mass displacement 5 CM of the IVC over a certain 
period of M MD steps 


5 
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Z ( r CM ( t m)~ r CM (0)) • 
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In the above expression, r CM (t m ) 


~N 


is the center of mass of the IVC containing N atoms at 

1=1 


positions r ^ at time t m of the m th MD step. 


In Eq. (1), the mass center displacement ?>cm * s calculated relative to the initial zero-displacement position of 
the IVC. This initial position is defined as the location of the mass center of the IVC when the MD system is 
equilibrated at zero applied load. 

The equilibration at zero applied load cannot be done using a MD domain with a free surface because the 
surface tension would produce a pressure, p s , on the surface resulting in erroneous zero-displacement positions for 
the IVC mass centers. To avoid surface tension effects, the zero displacements are estimated by using a MD system 
of a rectangular shape that is equilibrated under periodic boundary conditions in all three dimensions at zero 
pressure. The details of the equilibration procedure will be given in a separate publication 14 . 

Both the reaction forces calculated at the interface FEM nodes in response to the provided displacements and 
the MD interatomic forces are applied uniformly on the atoms of the IVCs as external forces. The application of the 
FEM reaction forces to the IVCs can generate undesired phonons or resonant elastic oscillations in the dynamic MD 
domain. These oscillations must be damped in order to achieve equilibrium with the static FEM region. A number 
of different damping schemes have been addressed in the literature. Holian and Ravelo 15 , and more recently, 
Schafer et al. 16 , found that applying viscous damping to the atoms in a region surrounding the center of the MD 
system can effectively absorb the intense phonon waves coming from a propagating crack. In this scheme, a friction 
force, fj , is given by 


b = -yy (2) 

and is applied to the atoms of the damped region in proportion to the atom’s velocity v and an appropriately chosen 
viscous coefficient x >5 ,16 . The method is efficient and simple to implement. Its drawback is that it effectively 
decreases the local temperature in the damping region resulting in undesirable strain gradients because of the 
thermal expansion. 

To avoid thermal effects, the damping used in the present method is based on a modified form of Equation 
(2), where instead of being proportional to each individual atom’s velocity, fj is set proportional to the group 
velocity of the mass center of a certain volume of N atoms. The frictional force, fj , and the group velocity of the 
mass center, v cm , are given by 


11 = -XV c 


1 N 
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By controlling the size of this volume, one can damp phonons of wavelengths larger than the volume size while 
leaving the shorter wavelengths associated with random thermal fluctuations unaffected as their contribution to the 
group velocity averages to zero. Phonons introduced by the FEM mesh cannot have wavelengths smaller than the 
distance / between the interface FEM nodes (Figure 2). Thus, it is convenient to choose the IVCs at the interface 
MD region and corresponding surface volume cells (SVCs, see Figure 2) at the surface MD region (to be discussed 
in Section ILC) as the volumes in which damping is applied. For the models used in the present work (to be 
discussed in Section III), effective viscous wave damping was achieved by setting % = 3 eVps/nm 2 . For this/, and 
for an IVC and SVC thickness of 1 and s, respectively, on the order of 3 nm, the effective temperature in the 
damping volumes decreased only by 1% compared to the bulk temperature. A discussion of criteria for the selection 
of the value of the viscous coefficient, j, can be found in Reference 16. In those instances where viscous wave 
damping is not adequate (e.g. collisional dynamics), the more precise non-reflective boundary condition techniques 
discussed in References 17-19 can be implemented. 

C. Surface MD region 

The proper embedding of a localized MD domain within a coupled FEM domain should generate a force 
environment within the inner MD region that closely approximates the internal forces of an infinite MD domain 
subjected to externally applied far-field loading. In ESCM, external forces are transferred from the surrounding 
FEM domain into the MD region as nodal reaction forces distributed to the atoms of the IVCs in the MD interface 
region. In order for the MD domain to deform freely in response to these reaction forces, it is simulated under free 
surface boundary conditions. The existence of a free surface introduces several undesirable effects in the MD 
system. First, it creates surface tension forces 
that must be removed to avoid distorting the 
MD response. Second, because atoms that 
should lie within the cutoff radius of the surface 
atom’s potential are missing, the surface atoms 
are less strongly bonded to the material than the 
bulk atoms. Under sufficiently high reaction 
forces, these atoms may be separated from the 
surface layer causing a surface degradation of 
the MD domain. To mitigate these free surface 
effects, a surface MD region - as shown in 
Figure 3 - is introduced. This additional region 
serves several necessary functions. First, 
positioning the outer free surface away from the 
interface MD region stabilizes the atomic 
structure of the interface MD region against 
strong FEM reaction forces produced during 
potentially large deformations. Second, the 
surface MD region is conveniently used as a 
volume where an additional external force is 
applied to counteract the surface tension force 
and to eliminate the surface tension effects on 
the remainder of the MD domain. 

However, this surface region introduces an undesirable fictitious stiffness outside the perimeter of the applied 
FEM nodal forces which elastically constrains the deformation of the inner MD domain. The combined effect of the 
surface tension and the stiffness of the surface MD region may be defined as a resultant force, f s , which acts at the 
boundary with the interface MD region and is given by the sum of two components expressed as 

/,=! + * (4) 

In this equation, £, is the elastic reaction of the surface MD region under deformation, and x is the force that results 
from the surface tension (Figure 3). 
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A way to mitigate both the surface tension and the elastic response of the surface MD region is to estimate and 
compensate the force / s . In the ideal case, when /" is fully compensated, the surface MD region acts as if it possesses 
zero stiffness and no surface tension, thereby mitigating spurious influences on the MD system. To follow the 
variations of/ s along the perimeter of the MD domain, the surface region is divided into a number of SVCs. The / s 
is calculated individually for each SVC. The force compensation is performed by calculating a counterforce, f c , 
along the IVC/SVC interface and applying this force over the atoms of each SVC in a similar way as the reaction 
nodal forces are applied to the IVCs of the interface MD region. Although the SVCs do not have to correlate in 
position and size with the IVCs, it may be computationally convenient to form SVCs as extensions of the IVCs into 
the surface region. The details of the force compensation method used are fully described in a separate publication 14 . 

D. MD-FEM coupling 

The MD-FEM coupling in the ESCM is achieved through an iterative equilibration scheme between the MD 
domain and the FEM region. In this scheme, the MD domain is simulated under constant traction boundary 
conditions, provided by the FEM as interface reaction forces Rp After a selected period of MD simulation, during 
which the MD domain responds dynamically to Rj, the displacements at the MD-FEM interface are calculated as 
statistical averages over the atomic positions and are averaged over time. 

These average displacements are imposed as displacement boundary conditions <5/ for the FEM region. The 
FEM problem is then solved to recover new interface reaction forces Ri, resulting from the interface displacements 
Si and the imposed far-field displacement boundary conditions Sp, which define the total loading conditions of the 
coupled system. The new interface reaction forces Ri are distributed to the interface atoms in the IVCs of the 
interface MD zone, defining new constant traction boundary conditions. The next period of MD simulation then 
takes place. Updated interface displacements are calculated and, again, provided as displacement boundary 
conditions to the FEM. This MD-FEM iteration cycle repeats during which a stable equilibrium of both 
displacements and forces between the atomistic and continuum material fields at the interface is established. The 
period of MD simulation prior to FEM update of Ri is selected by a determination of the convergence rate to a stable 
equilibrium state. 

The basic equations of the ESCM coupling methodology are now developed. The elastic properties of the 
material in the FEM region are described through a set of stiffness matrices \K u y\. with a, ft = V, F, 1 indicating: V — 
variables within the interior of the FEM region; F - far-field variables; and / - interface variables. Using these 
definitions, the static continuum equations of state at the /7 th FEM update are given by 


'[Kyy] 

[ k vf \ 

[Kyi] 




[ k fv\ 

[ k ff ] 

[ k fi ] 

' { 5 f(0} 

> = < 
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To ensure that the FEM region has the same elastic properties as the MD domain, the stiffness terms in the 
[ K u y \ submatrices are calculated from the anisotropic elastic constants derived from the MD interatomic potential at 
the temperature of the simulation 20 . The FEM model is subjected to two types of displacement boundary conditions: 

(i) the far-field displacements {$f} , which define the load over the entire coupled MD-FEM system; and (ii) the 

interface displacements {§/} = (§/,§/,.. §/,..) , which represent the deformation response of the MD domain at the 
1 st , 2 nd , ... , and k th IVC. 

The solution for the unknown displacements in the interior of the FEM domain, {5f/} , is given by 


{*V(0} - \ K vv\ ! (K } - \ k vf }{S f ( c n )) + [ K vi r (hi )} ) > 


(6) 
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which allows the calculation of the interface reaction forces, {/?/(?„)} = (Rj,Rj ,.. Rj ,..) of the 1 st , 2 nd , ... , and k th 
IVC to be obtained from 


{ R i {h i )} - 1 , K iv ]{§ v n )} + [ k if ] 
together with the far-field forces of constraint 

FV ] MOM * ff\ 


n )\ + \ K II ]{§/ ( t n )} > 


{5F( ? n)}+[^F/]{5/(b;)} • 


(7a) 


(7b) 


The dynamics of an atom i of mass |i (l1 at position r (l) in the embedded MD domain is described by Newton’s 
equations of motion 


|a (0^(0 = y(0. ;'e (inner MD region) 

p (0^; 0)_y(0 +r^/Nj -j\ k cm \ is(lVC)/. interface shell 
H U) r (,) = f (,> +fc/ N s ;e(SVC ) k surface shell 


( 8 ) 


where % and V cm are defined in Equation (3). 

Convergence of the interface forces {/?/} and displacements { 8 \ | can be achieved only if the MD system can 
reach a force balance with the exterior FEM system subject to the far-field displacement constraints. At equilibrium, 

first, x k cm = 0 , and second, the change of the total momentum Ap k of any k' h IVC due to the FEM reaction force for 
the period of the MD run of M time steps is also zero 

M N I 

A Pk = 'L Z ft {l) r (,) (t n )At = 0 . (9) 

n = 1 i=l 


Performing the same double summation to the second equation in (8) results in 


, M N i 
M m=l /'=! 





(10) 


which means that, at equilibrium, the FEM reaction force becomes equal and opposite to the average MD atomic 
force in the corresponding /c th IVC. Equation (10) presents the establishment of static force equilibrium between the 
MD domain and the FEM region. 


III. Test Simulations of the ESCM 


A. The simulation models 

A series of tests is presented in this section to examine the performance of the ESCM. The model geometry 
used in all the tests is shown in Figure 4. This model consists of an inner circular MD domain of diameter c/vm 
which is embedded into a larger exterior square FEM domain of elastic material with a side dimension of d vv = 
20c/ M |). The specifics of the MD and the FEM models and simulations will be discussed next, followed by the 
results and analyses of the test simulations. 
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Figure 4. Coupled MD-FEM model. 


1. The MD model 

The material for the simulation models is chosen to be a perfect crystal of aluminum. The atomic properties of 
aluminum are represented by the embedded atom model (EAM) potential of Mishin et al. , which is fitted to give 
the correct zero-temperature lattice constant, a = 0.405 nm, elastic constants, cohesive energy, vacancy formation 
energy, etc. A 2000 point tabulation of the potential function in the interval of interatomic distances from r min = 
0.25u o (0.1 nm) to the cut-off radius, r c = 1.55a 0 (0.628 nm) was used. 

In this study, a second order polynomial interpolation for the tabulated EAM potential was applied, which, 
combined with a high precision fifth order Gear predictor-corrector integration scheme for calculating the atomic 
trajectories 22 , gave errors in estimating the stress in the system under uniaxial strain in the range of 0.1 - 0.5% that 
were less than 0.5 MPa. This is equivalent to detecting a change in strain of the order of 0.007%. The high accuracy 
in strain detection allowed accurate measurements of the MD-FEM coupled model at relatively small strains of 0.1 
to 0.5%. 


The first test simulation used a MD system in the form of a circular plate of monocrystalline aluminum with its 
main crystallographic axes [10 0], [0 1 0] and [0 0 1] oriented along the x-, y- and z- directions, respectively 
(Figure 4). The system was simulated under periodic boundary conditions along the z- direction and free surface 
boundary conditions along its perimeter. These boundary conditions allow the MD domain to deform in an 
unconstrained manner in the x-y plane under the external reaction forces from the FEM region while maintaining 
constant zero pressure along the z-direction using the Parrinello-Rahman constant-pressure simulation technique 23 . 
Constant temperature is maintained by applying the Nose-Hoover thermostat 24 . The tests were performed at near 
zero temperature (T= 10 K) to minimize the thermal noise, and at room temperature (T = 300 K) to demonstrate the 
ESCM for more practically relevant situations in which thermal effects are important. 

The thickness of the plate along z is equal to the thickness, h = 5 a 0 ~ 2.0 nm, of the coupled MD-FEM system as 
shown in Figure 4. Though very thin, the MD domain mechanically behaves as an infinitely thick plate due to the 
applied periodic boundary conditions along the z- direction. The interference of the atoms with their periodic 
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images in the z-direction is prevented by the fact that h is more than three times larger than the range of the 
interatomic potential r c thus preserving the local three-dimensional characteristics of the system. 

Although the ESCM approach can be applied to an arbitrary interface geometiy, the circular shape of the MD 
domain was chosen for several reasons. First, the circular shape is the equilibrium shape of a free body with respect 
to the surface tension. As already discussed in Section II, in the ESCM approach the MD system represents a free 
body placed in the elastic force field generated by the continuum FEM domain. The MD-FEM coupling is more 
effective if the free body is in equilibrium with itself before the coupling is imposed, which is not guaranteed for an 
MD domain of general shape. Starting with a nonequilibrium shape of the MD region may cause the applied 
counter-forces to be significant compared to the interatomic forces and, thus, may significantly affect the atomic 
state. Second, the MD-FEM interface is uniform everywhere around the MD system, which significantly simplifies 
the analyses of its behavior and the surface tension compensation procedure as described in Section II. C. Third, the 
curvature of a circle (or sphere) undergoes minimal changes during deformation, compared to any other shape 
because the radius of curvature is constant and equal to the maximum possible value everywhere on the surface. 
Thus, the change of this radius under strain will be minimal. Preserving a constant curvature of the free surface 
during deformation causes the surface tension to remain constant during the simulation, which additionally 
simplifies the surface tension compensation procedure. 

2. The FEM model 

The elastic continuum region was modeled using 8-node hybrid-stress hexahedral finite elements 25 ' 26 that have a 
reduced sensitivity to mesh distortion compared to standard displacement-based elements. The constitutive relations 
were selected as linear elastic to avoid the uncertainty of introducing nonlinear terms in the material constitutive 
matrix, which are not trivial to estimate from the interatomic potential. The elastic constants in the material 
constitutive matrix were based on the MD-estimated elastic constants of pure aluminum at T = 10 K. The values 
were averaged for uniaxial stresses from 100 to 500 MPa as: Cn = 112.7 GPa, C [2 = 59.4 GPa and C 44 = 30.6 GPa. 
These values are about 3% different from the static, zero Kelvin, elastic constants that have been reported for this 
potential 21 . 

The continuum finite element model, as shown in Figure 4, contains an open inner region of diameter d M D , in 
which to embed the atomistic domain. Along its perimeter, 80 nodes at z = +/?/ 2 and 80 nodes at z = -h/2 were 
placed to form 160 FEM interface nodes to communicate with the embedded MD domain. 

The dimensions of the FEM mesh, d\v, d\ m and h (see Figure 4) were defined initially through their proportions 
d EE '■ duet : h = 20 : 1 : 1 . Then, the dimensions <7md and h were rescaled to the equivalent dimensions of the MD 
system after MD equilibration at zero stress and constant temperature, and after the width of the MD-FEM interface, 
/, and the surface MD region, s, had been determined. Finally, the outer dimension of the FEM system d EE was 
defined to preserve the ratio d EE : d MD = 20 : 1 . The outlined procedure allows a single FEM mesh to be applied to 
various MD systems of different structures and crystallographic orientations through geometric rescaling. This 
universality of the FEM mesh in the ESCM is a convenient advantage compared to the direct coupling methods 1 " 11 , 
in which the FEM mesh has to correlate to the atomic positions of the interface atoms and is specific for each model. 

B. Test results and analyses 

Two test cases are considered to interrogate the behavior of the ESCM. First, the stress-strain continuity 
between MD and FEM domains is validated using the model of an elastically deformed plate with a circular hole. 
Second, simulating the propagation of an edge crack through the FEM domain into the MD domain is performed to 
gauge the overall representational capability of the coupled MD-FEM model. 

1. Test for stress-strain continuity over the MD and the FEM regions: elastically deformed plate with a 
circular hole 

To test the efficiency of the ESCM in transmitting the stress-strain field of a large elastic continuum region to a 
much smaller MD atomistic region, the classic example of a plate with a circular hole subjected to uniaxial loading 
was used. In addition to having a simple exact elasticity solution, this model is particularly convenient for two 
reasons. First, it provides large stress variations (from zero to 3cr 0 ) around the hole, which can be used to test the 
continuity of the stress field at the MD-FEM interface in the case of large stress gradients. Secondly, keeping the 
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peak stress, 3cr 0 , well below the theoretical elastic limit of the material (recently estimated for the applied 
interatomic potential to be between 3 and 5 GPa 30 ) prevents the occurrence of any plastic mechanisms in the MD 
domain, which is not addressed in the present study, and a static elastic equilibrium can be achieved everywhere in 
the system. 

To develop the model, a hole of radius 20 nm is removed from a FE model of dimension 800 nm by 800 nm and 
having properties matching the MD elastic properties for a perfect crystal of aluminum 21 . The elastic anisotropy of 
the plate follows the crystallographic orientation of the atomistic monocrystalline lattice in the MD domain given in 
Section III. A. 1 . The aluminum material from 20 to 40 nm from the center of the hole is simulated atomistically by 
a circular MD region with diameter c/ MD = 80 nm (200«„). The MD region is embedded in a continuum elastic square 
plate of size < f FE =1.6 pm (Figure 4). The square plate is deformed at 0.5% uniaxial strain along the x-direction 
through fixed displacement boundary conditions imposed on the outer sides of the FEM system 800 nm away from 
the hole in the MD domain. At this strain, the far-field stress, estimated as in Figure 5, is a 0 = 358.5 MPa. The 
stress components in the z-direction are kept at zero by controlling the thickness h of the plate using the constant 
stress Parrinello - Rahmann technique applied along the periodic z-direction in the MD simulation. 

Starting from an undeformed MD domain, |5 j } = {o}, the coupled MD-FEM iterative simulation is performed, 
as described in Section II. D, until equilibrium is established between the MD domain and the FEM region ( {5/} and 
{Rj} converge to constant values). The stress 

field for (T cv , < jyy , and a xy stress components of fem,-md fem 

the coupled MD-FEM system is calculated and 
compared with a full three-dimensional 
anisotropic continuum FEM solution, 
performed using the ABAQUS software 
package 28 (Figure 5). 


The stress profiles for a xx and a yy along the 
x- and y-axes calculated through the center of 
the hole are shown in Figures 6(a) and 6(b), 
respectively. The continuity of the stress 
profiles is well preserved at the MD-FEM 
interface apart from fluctuations that result 
from chaotic atomic thermal vibrations. The 
magnitude of the thermal fluctuations can be 
seen in the simulation data of a continuous 
plate corresponding to the MD-FEM model 
shown in Figure 5 but without an open hole and 
subjected to the same far-field strain (cr rx (p i a , e) 
and <7,y ( p i ate ) in Figures 6(a) and 6(b)). In 
Figures 6(a) and 6(b), the symbols represent 
the ESCM simulation data, while the lull dark 
lines represent the fully continuum ABAQUS 
simulation data without surface tension. The 
stress profiles of the coupled MD-FEM model 
follow closely the stress profiles of the fully 
continuum FEM analysis of the same model. 

The largest discrepancy is observed at the 
edges of the circular hole. One reason for this 
discrepancy is the surface tension that exerts a 
tensile force towards the center of the hole in 
the MD simulation. Another reason may be 
because the virial stress gives poor 

convergence and erroneous estimates near free 
surfaces 29 Figure 5. Comparison of stress calculation in MD-FEM and 

FEM models. 
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Figure 6. Comparison of stresses between MD and FEM models along cuts in the x and y 
directions through the inner open hole. 


In an attempt to account for the surface tension of the hole in the continuum model, an external pressure of p s = 
45 MPa (estimated as p s = yfr with y= 0.9 J/m 2 , being the energy of the free surface 21 of the hole, and r = 20 nm, 
being the hole radius) was applied normal to the surface of the hole in the FEM model. Unlike the continuum 
simulation without surface tension, the model with surface tension correctly represents the normal component of 
stress at the hole’s surface, which becomes equal to p s and is plotted using grey lines in Figures 6(a) and 6(b) 
(following <r xx along the X- pro file in Figure 6(a), and a rr along the y-pro file in Figure 6(b) when approaching the 
hole radius at ±20 nm in both plots). The corresponding MD profiles are in good agreement with the continuum 
prediction with surface tension. In contrast, the tangential component of the stress (<r n along the x-profile in Figure 
6(a), and a xx along the y-profile in Figure 6(b)) in the MD system deviates substantially from the continuum 
prediction when approaching the surface of the hole. As shown in Figure 6(b), the continuum solution closely 
approaches the theoretical value of 3 cr 0 that is calculated for an isotropic material (note that the continuum model in 
this simulation is anisotropic and some deviation from the isotropic solution is to be expected). As mentioned in the 
previous paragraph, one reason for this discrepancy may be the virial definition of stress, but more likely it is that 
the continuum model does not account correctly for the nature of the surface tension, which results from the 
occurrence of missing bonds between the atoms at the free surface. The effect of missing bonds may result in much 
stronger local tension between these atoms at the edge than can be considered in the continuum simulation. 

2. Example of an edge crack simulation along a grain boundary in aluminum 

A possible application of the ESCM technique described in this paper is the simulation of atomistic processes 
where satisfying the correct elastic boundary conditions is of a crucial importance. A typical example is the tip of an 

edge crack, where the idealized elastic stress field decreases as 1/ sfr and extends to distances r from the tip far 
larger than a MD domain can computationally simulate. A coupled FEM-MD model for this example is presented in 
Figure 7, where the MD domain is used to atomistically simulate the plastic zone at the crack tip, and the FEM 
region is used to provide the continuum elastic boundary conditions of a far-field tensile strain of £ yy = 2% applied 
along the v-direction. To check the validity of the technique at higher temperatures, the simulation was carried out 
at room temperature (T= 300 K). 

The dimensions of the system are: h = 2.9 nm, duo = 45 nm, and <7 fe = 900 nm (see Figure 4). The MD domain 
represents a bicrystal with a high-angle grain boundary formed between the two crystals, along which the edge crack 

propagates. In the imposed coordinate system of the model, the orientation of one of the crystals is: (x:[ 7 7 10 ]; 
y:[5 5 7 ];z:[l 1 0 ]), and the orientation of the other crystal is: (x:[ 7 7 10];y:[5 5 7 ]; z: [ 1 1 0 ]). 
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In this way, both crystals are mirror images of each other relative to the crystallographic plane {5 5 7 }, which 


What is essential in the example, provided 
in Figure 7, is that the developed ESCM 
approach preserves the continuity of the stress 
- strain field at the MD-FEM interface even 
for a dynamic problem such as crack 
propagation simulated at room temperature. 

The depicted profile shows a good continuity 
of stresses across the boundary when the crack 
speed is slow compared to the elastic response 
time of the system. In the example shown in 
Figure 7, the crack propagation speed was 
around 100 m/s or of the order of 1/30 the Figure 7. Simulation of an edge crack using ESCM. 

speed of sound. 



becomes the plane of the grain boundary (GB) between them. The GB thus formed is a <1 1 0> £99 symmetric tilt 
GB. Crack propagation along this GB has been extensively studied by MD simulations 29 ' 31 at a cryogenic 
temperature of T = 100 K. In these studies, it was found that while in one direction the crack becomes blunted by 
deformation twinning, it propagates in a brittle-like manner in the opposite direction with very little dislocation 
emission. The latter direction has been chosen as the propagation direction for the edge crack of the simulation in 
Figure 7. The simulation using the ESCM approach performed at room temperature (Figure 7) showed higher 
dislocation plasticity than at T = 100K 30 " 32 . The problem remains how to transfer this plasticity to the FEM 
continuum. Some work in this aspect has been 
started by Shilkrot et al. s where they have 
developed a method to transmit a dislocation 
from the atomistic region to the continuum 
region and to follow its propagation using 
dislocation dynamics simulation in 
combination with the FEM. 


Continuum field 


45 nm 


IV. Concluding Remarks 

A new statistical approach to coupling MD with FEM simulations denoted embedded statistical coupling 
method (ESCM) is developed. The approach is based on solving the boundary value problem through an iterative 
procedure for both MD and FEM systems at their common interface. The two systems are simulated independently 
and they communicate only through their boundary conditions. The FEM system is loaded along the MD-FEM 
interface by nodal displacement boundary conditions obtained as statistical averages of the atomic positions in the 
MD system at the mass centers of associated interface volume cells (IVCs). The MD system is simulated under 
constant traction boundary conditions obtained from the FEM system as reaction forces to the MD displacements at 
the interface. This approach allows the continuity at the MD-FEM interface to be achieved at different scales 
inherent to both systems without the need of redefining continuum quantities at the atomic scale and atomic 
quantities at the continuum scale, thus avoiding a theoretically controversial concept. 

The ESCM shows excellent convergence for systems simulated under static elastic equilibrium, and preserves 
stress continuity at the MD-FEM interface for systems exhibiting relatively slow dynamics. Additional study needs 
to be performed to verify if the implementation of a dynamic FEM simulation can improve the dynamics of the 
system away from equilibrium. 

Compared to direct coupling methods 1 " 11 , ESCM does not represent the MD-FEM interface down to atomistic 
dimensions but uses statistical mapping of atomistic behavior to continuum FEM representation. ESCM presents a 
simple and flexible technique in providing elastic boundary conditions for an MD model, thus solving to a large 
extent the finite-size artifacts inherent in an atomistic simulation. By using statistical averaging to reduce the 
degrees of freedom when passing information from the atomistic system to the continuum model, the FEM mesh can 
be independent of the atomic structure of the MD system, and the same FEM mesh can be used with MD models of 
different size and microstructure without the necessity of remeshing during the simulation. This may prove a 
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helpful advantage when plasticity is considered. There is no limitation on the temperature of the atomistic system as 
long as the thermal corrections to the elastic properties of the continuum system are considered. 

From a computational point of view, there are several differences between the ESCM and DC approaches. First, 
because the FEM mesh is not refined down to the atomic scale at the interface, a significant computational 
advantage can be obtained when compared to existing DC methods because the size of the finite element model can 
be decreased by orders of magnitude. Thus, computational resources can be focused on the atomistic calculations 
yielding the possibility of embedding large MD domains on the order of millions or billions of atoms within the 
overall material model. In general, an increase in the size of the MD domain associates more atoms with each FEM 
node in the interface region and improves the continuity at the MD-FEM interface due to improved statistical 
averages. The integration of detailed atomistic resolution and relatively coarse refinement of the FEM continuum at 
the interface enables the construction of large models that may benefit engineering applications in the area of 
nanomechanics such as the study of nano-electromechanical systems (NEMS), integrated circuits, and thin film 
structural configurations. 

Second, because the mesh does not have to correlate with the atomic structure at the interface, it may be used 
with different MD domains of different atomic microstructure. This independence of the microstructure could be 
very useful in considering plastic processes at the interface, as no remeshing is needed if atoms experience plastic 
flow. Also, from a practical standpoint, universal meshes can be prebuilt with a fixed interior region used to embed 
the MD domain. Thus, various MD systems used to study different atomistic processes within the same overall 
material model can be repeatedly used without having to construct a new FEM mesh for every different atomic 
configuration. 

Additionally, the method is relatively easy to implement. Any conventional FEM code, including commercial 
packages, can be used to solve the FEM part of the model. Small modifications to perform the constant traction MD 
simulation (outlined in Section II. A) are necessary to an otherwise conventional MD code. Thus, the method may 
be suitable for use in nanoengineering applications where small regions of special interest can be simulated at the 
atomic level while coupled to a larger material system that can be represented as a continuum. 
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